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FORMULATION AND ERROR 
ANALYSIS FOR A GENERALIZED 
IMAGE POINT 
CORRESPONDENCE 
ALGORITHM* 

SUNIL FOTEDAR 

Advanced Automation Technology Group , McDonnell Douglas Corporation , Space 
Station Division, Houston, Texas 77062 

RUI J. P deFIGUEIREDO 

Department of Electrical and Computer Engineering , University of California , Irvine , 
California 92717 

and 

KUMAR KRIS HEN 

Johnson Space Center, NASA, Houston , Texas 77058 


A Generalized Image Point Correspondence (GIPC) algorithm, which enables 
the determination of 3-D motion parameters of an object in a configuration 
5 where both the object and the camera are moving, is discussed. A detailed error 
analysis of this algorithm has been carried out. Furthermore, the algorithm was 
tested on both simulated and video-acquired data, and its accuracy was 
determined. — _ _ 


I. Introduction 

Motion analysis, based on robotic vision, has widely been discussed in the 
literature [1-5] in developing the Image Point Correspondence (IPC) algo- 
rithm. Methods used are Two-view motion analysis or monocular vision , stereo 
or binocular vision , and stereo motion . However, the IPC algorithm has not 
been applied to the more general problem of motion analysis involving a 
situation where both the object and the camera are moving [4]. Industrial 

•This research was supported by NASA Contract NAS 9-17145 and NASA/JSC (RICIS) 
Grant NCC 9-16. Thanks are due to Mr. Olin Graham. Special thanks to Mrs. Lovely K. 
Fotedar for technical assistance. 
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and space robots face this situation in locating and tracking various objects/ 
scenes when both the camera/video system and tlje object move asynchron- 
ously. In the GIPC, a generalization of the IPC algorithm, this problem of 
motion analysis is discussed. The three methods of motion analysis men- 
tioned before become special cases of the GIPC presented in this paper. 

The accuracy of the IPC/GIPC algorithms depends on the availability and 
errors in the measured input parameters. Input errors also include the 
deformation of the object over time. The error propagation for various stages 
of the IPC will, therefore, be described in terms of the error bounds. 


II. The Generalized Image Point Correspondence 
Algorithm 


2.1. The Algorithm 


PAGE is 

OF PO OR QUAU1 i 


The general case of motion analysis is illustrated in Fig. 1. For simplicity in 
presentation, the equations that track a single point P on a moving object 
viewed by a moving camera are considered. F f and Fj are the two frames with 
which the camera coordinate system coincides at two different instants of 
time tj and tj ( > r f ), respectively. Point P moves from one position P { to 
another position Pj due to the rigid-body motion of the object. We assume 
(R h T f ) and (R p T ; ) to be the transformation parameters (rotation and 
translation) that link the frames F* and Fj respectively with the standard 
frame S. Also, let (R ijy T lV ) be the transformation parameters that link the 
frame F f with the frame Fj. The object moves with the unknown motion 
parameters (R, T). The image plane is assumed to be at the focal point of the 
camera with its X - and T-axes parallel to those of the camera coordinate 
system, where Z-axis is the line of sight. . 

The desired relationship between the coordinates of the initial and the final 
positions of point P ( P f and P jy respectively) recorded by the camera, with 
respect to the frames F, and F jy respectively, is given by the equation [4] 

Pjj = R'u Pa + T„ f (2.1a) 

where p a/I = (x afl , y afiJ r a ^) T is the vector of 3-D coordinates of P fi relative to F m 
at instant t p (a, /? = L j\ and 


R' u = Rjj Ri RRJ = Rj R Rj and T u = - RjR RjT { + R,T + T,. 

(21b, c) 


Equation (2.1a) gives the expression for the generalized version of the 
motion-analysis equation. Clearly, it does not matter whether the object or the 
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OBJECT UNDERGOING RIGID-BODY MOTION 



camera is moved first. R and T, the desired parameters to be estimated, are 
defined as 


>11 

r 12 

r n 


v 

'21 

r 2 z 

r 23 

and T = 

*2 

r 31 

r a 

r 3 J 


h 


(lld,e) 


where r, f (a, p = 1, 2, 3) are the rotational elements and r, (a = 1, 2, 3) the 
translations along x-, y-, and z-axes, respectively. 
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Special Cases. If F, is standard frame S, the motion equation can %e 
written as 

R'ij =R]jR = RjR and = RjY + (2.2a,b) 

The I PC algorithm can be used to estimate the motion parameters 
(Ry, T; v ) and hence (R, T) of the moving object, assuming R, v and are 
known. 


2.1.1. Monocular Vision 

In the case of the two-view motion-analysis equation, the location of the 
camera taking the pictures of the moving object, is fixed. In that case, 

Ri = Rj = R u = / and = T, = O. (2.2c,d) 

The generalized motion equation, using Eqs. (2.2a,b) and (2.2c, d), reduces to 

P jj = RPii + T 

or to the more familiar two-view motion equation 

p' — Rp + T 

as the frames F f and F j coincide with the frame S, such that p a = p { = p and 
p jj = p j = p', where p t = (x ; , y h z,) T and p j = (x ; , y jy Zj) T are 3-D coordinates 
of P ( and Pj relative to S, respectively. 


2.1.2. Stereo Vision/Stereo Motion 

For a stereo vision/stereo motion case, the object is assumed to be 
stationary. In that case, 

R = /; R u = Rj and T = O (2.2e4) 

and the generalized motion-analysis equation, using Eqs. (2.2a,b) and (2.2eJ), 
reduces to 

P' = RjP + Tj 

These cases of motion analysis have been found to be equivalent [1]. The 
motion of point P t to point Pj with respect to a fixed frame Fj is similar to the 
motion of frame Fj to frame F ( - with respect to a fixed point P. 


owgimi ;.\^t is 
OF POOR QUALITY 


III. Error Analysis 


Since the IPC/GIPC algorithms can be implemented using a sequence of 
object images, any error in the input data and sensor parameters becomes a 


\ 
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Fig. 2. Error analysis for various stages of the IPC algorithm. 


source of inaccuracy in the output data. The input data are the set of feature 
coordinates of the object, and the output data are the 3-D motion param- 
eters. The sources of perturbation errors have been identified in [6, 7], 

In this section, we analyze the effect of changes in the input parameters on 
the output parameters for the three different methods for the IPC/GIPC 
algorithm [1-4]. The two representations of the rotation matrix will also be 
considered [2, 4, 8]. Since various steps are involved in this algorithm, the 
propagation of the error will be studied and the error bounds for each stage of 
the algorithm (Fig. 2) [7-9], 

3.1. Error in Essential Elements 

The sensitivity of the essential elements (elements of a matrix Q) to the 
error in input set of 2-D image coordinates of the features, common to the 
three methods, will be studied in this section. By definition, matrix Q is 
represented in terms of rotational and translational elements as [2, 4] 


Q = 

9n 
9 21 

<?12 <7l3 

<?22 <?23 

= ± 

*3 r 21 ~ *2 r 31 
*l r 31 “ h r n 

*3 r 21 *" *2 r 23 *3 r 2l f 2 r 33 

^l r 32 “ *3 r I2 *l r 33”*3 r X3 


9m 

<?32 <?33 


h r n — *l r 2l 

f 2 r 12 “ *I r 22 *2 r 13 "" *l r 23 


(3.1a) 
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In our case, the equation solved is 

- 1 )* 

and the true equation to be solved is 

= -i. ...,-D T , 

where N = N + 5N; Q = Q + 5Q; 

N = (N lt N 2 ,N 3 ) T ; 

n. = (a,a;, a-;, a;, a., y„, r„ r„ a., a.) 

> 8 ); 

and each element of Q is divided by q 33 . Let 

5N' = (5a xl , 5 a x2 , Sa x3 da x9 ), 

where 

Sa x , = <5 AiA. + <5A.A; ; Sa x2 = 5X' x Y a + 5 Y x X' x , 

So . 4 = SY' x X x + SX xl Y' x ; 6a xi = + SY x r x , 


(a = t,2...,n;n 


5a. 3 = <5A; : 
&*«« = <5 a;,; 


&J a7 = <5A,, and <5a >8 = £A a . 

The products of errors are neglected and (X xfx ) are 2-D image coordin- 
ates for the ath data point. It follows from [9] that 


N SQ = — SN Q = — 


■*r 

* 


* 





Sa 1 1 ^ t Sa | g 

* * * * 

• * * * 

* * <5<*o8 



9n ' 


* 


* 


«J2. 


. (3.1b) 


Thus, if the inherent errors 5a xf were known, the corresponding solution 
errors Sq lf would be obtained by solving Eq. (3.1b). The degree of accuracy is 
consistent with the assumption of neglecting product of errors. Based on 
dimensions of N, we have the following cases: 


Case /: If N is square (i.e., n = 8) and non-singular, 

\*1*\£E* (« = 1 , 2 ,..., 8 ), 

where 


(3.1c) 


E. = «« Z Z (a, P = 1, 2, 3; a and 0 * 3 => q 33 not included) (3.1d) 

or = 1 0- 1 


for -c a <. Sa xf < i x (<x,p = 1,2,. . 8). Therefore, 

8 

\Sq^\ < E a Z I A 3(a _ „* # .,| (for a, 0 = 1, 2, 3; a and /? # 3), 


(3.1c) 
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where (a, fi = 1, 2,..., 8) are the elements of the (Ilh row of N - I. Thus, if 
the elements of AT 1 are calculated, approximate upper b6unds on the effects 
of inherent errors can be obtained from Eq. (3.1e). Since* they were derived 
under the assumption that the products of errors are negligible relative to E a , 
they are not strictly upper bounds. However, they are acceptable as clow 
approximations to the true upper bounds. 

Case II. If N is not a square matrix, we define a residual matrix by 

5N + = + - N + , 

where N + and f)+ are the pseudo-inverses of N and f} respectively Then 
from [10] 

II-SATII < mt\\ + U 5 N 2 || + \\6n;\\, (3.10 

where 

«^. + ! ^ ll^ll • II W + II • II N+L WSNIW Z |Mf| • and 

ll<5M 3 + || < \\SN\\N + f. 

3.2. Error in Rotational Elements 

The matrix Q is found in all the three methods in the same manner. But the 
rotation matrix R is computed from Q in three different ways. In this section, 
we shall treat these methods separately in order to investigate the propaga- 
tion of error due to error in essential elements. 


3.2.1. Error Analysis using the First Method 

The propagation of error in the rotational elements, when the first method 
for the IPC algorithm is used [2,4]. is investigated in this section. If the 
singular value decomposition (SVD) of Q, defined as 

Q = Q + sq = Gav t = (u + SUX A + MXK+ SV) T , 
is an approximation to the true value 


Q = U A V T , 


where SU , M, <5 Pare the error matrics, then the ath singular vector fi* of the 
perturbed matrix Q in terms of the ath singular vector u„ of the matrix Q can 
be approximated by the first-order Taylor Series Expansion and is given by 


du 

“■ = u ‘ + nT 

d %y 


1 4*v~**r 


(9#r-W (M.y= 1,2,3). 


(3.2a) 
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Therefore, 



du_ 


* 

IKI £ 

« 

8 %y 

1fiy = <lfiy 



(3.2b) 


where q fly , q fiy are the entries of the matrices Q and (J, respectively. The 
derivatives of the singular vectors (for k = l, 2, . . 9), assuming the singular 
vector V is known, are 




10 y = 4#y 




X [-r^-r (“I v,) + - « B « 0~k (3.2c) 

°yy ^yy) J 

7“l 


' - {o. 


1, if a + /i — 1 = k = 3(/J — 1) + y 
otherwise 


(3.2d) 


where <r„ (a = 1, 2, 3) are the singular values of Q. The advantage of using 
Taylor Series Expansion is that it provides the first-order derivative of the 
principal singular subspace with respect to the essential elements, and it can 
be combined with other derivatives via the chain rule to obtain first-order 
perturbations in essential parameters. In this approach, the singular vectors 
are considered to be vector-valued functions of essential elements. This 
approach, however, works only if we have distinct singular values. The 
expressions become ill-conditioned as singular values get close together. 

After finding error bounds for the singular vectors, we are in a position to 
find the same for the rotational elements, which are given by the equations 


I (K, + Mr I) (a, /J = 1, 2, 3), (3.2e) 

7=1 

where u afii o afi , and v afi are the entries of (/, A, and K respectively, and 

£ ucu = = £* = £ v for ~ £ u ^ ^ ^ ^ and ^ ^ 

e v has been assumed. 


3.2.2, Error Analysis using the “ Second Method" 

The rotational elements in terms of essential elements, using the second 
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method [1,4], are given as 

4 

r a,fi — ± [<?*,0 + 2OZ a + 1.0 + 2<7«+2,0 ~ 1.0<7« + 2,0 + 2) 4 

— *7a,0 + lO?at + l,0<?« + 2,0 + 1 ~~ <f* + 1.0+ l^a + 2,0) 

3 3 

— *«(?« + l.0 + l<7«+2,0+2 ~ <7*+1.0 + 2<7*+2,0+l)]A a Z Z q a 0> 

ar* l 0- I 


(3.3a) 

where a, /? = I, 2, 3, and are cyclic (for instance, for a = 2, /? = 3, 
^8 + 2,0 + 1 = 9u)« 1 * ^ 2 * * 3 ) arc translations along the x-, y-, and z-axes 

respectively. Therefore, the error bounds for the rotational elements are 


3 3 


i^i^, z z 

7 ~ 1 IC* 1 


<>%* 


(3.3b) 


f° r ~ £ q ^ ^ E ,i (y, k = 1, 2, 3). The partial derivatives in Eq. (3.3b) have 

not been computed in this paper. 


3.2.3. Error Analysis using the “ Third Method f" 

In this section, the error propagation using the third method for computa- 
tion of R from Q [3] is presented. Defining H = Q r and 

W. = h,xT = W.,f + W c2 j + W. 3 £, (3.4a) 

where h„ is the ath row of H, (i, j, £) are the unit vectors, and 

= (*lt+i.ttfi+2 ~ Qf+2,*tf+i )> (3.4b) 

where a, /? = 1, 2, 3, and are in cyclic order. Therefore, rotational elements are 
given by 

r *f Kfli+2,f+i^ 9i.rti^+2Xfl^r+j^+i - i.i+jU 

— (tfot.f+l^+l — + ~ 9«. 0 + 2*« + 2 )] 

+ C9«+i.#*«+ 2 — 9«+2.# t «+i]- (3-4c) 


Hence, the bounds are defined as 


l&Vl = e. 


3 


I 


d<lr« 


for — £, ^ Sq y „ <. e f (ot, /?, y, k = 1, 2, 3). (3.4d) 


The partial derivatives in Eq. (3.4d) are not computed in this paper. 
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3.3. Error in Translational Vector 

4 

The translational elements f a (a = 1,2, 3), in terms of essential elements, are 
given as [2] 

f, = /Z (-<& + <?2 +i.* + <72+2.#) (a is cyclic). (3.5a) 
V 0 = 1 


The error bounds for these elements are given as 


l*.i E Z 

0-1 7=1 


q Jl 

f. 


for — < <5, yi[ < e, (y, k = 1, 2, 3) because 


K 

d< fl>y 


+ ^ for /l# a 

I* 

— — for p = x 


(3.5b) 


(3.5c) 


3.4. Error in 3-D Motion Parameters 

In this section, the errors in the motion parameters due to errors in 
rotational elements, are studied separately for two representations of the 
rotation matrix R [2, 7, 8]. 


3,4.1 . Using the First Representation of the Rotation Matrix 


From the definitions of the directional cosines of an arbitrary axis, and the 
angle of rotation around this axis using the first representation of R [2, 7], the 
error bounds for these motion parameters are found to be 


3 \d v | 3 

|<5v a | < c r £ k- 1 and \8B\<.t r £ 

0.7=1 \ or fy\ a.0-1 

0* 7 


30 

dr.. 


(3.6a,b) 


for — e q ^ S rfiY ^ e r (a, /?, y = 1, 2, 3). General expressions for the partial 
derivatives are 


dVq 

5r 0 7 


( r o + 2.0+1 r o+ l,o + 2X r 70 — r ty) 

d 3 


, d l ~ ( f 70 ~ f 0>) 2 
~ d 5 


for /?, y = 1, 2, 3; /? # y and 
a is cyclic 

for p,y ^ a; P,y in cyclic order 
(3.6c) 
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and 


89 

dr* f 


± ( r *! ~ r e*) 

dj 4 - ~d 2 


I 

(for a = I, 2 , 3; /? = a + 1 ; a,*/? are cyclic). ( 3 . 6 d) 


3.4.2. Using the Second Representation of the Rotation Matrix 


The second representation of the rotation matrix R is used in this section 
[ 6 , 8 ]. The error bounds for the motion parameters in this case are: 


where 


|50|<; 


dO 

dr zs 


l^l 



m 

< 

dtp 
dr ,j 

\dr I3 | + 

8 <p 
dr 3 3 

1^331 


m 

< 

dij/ 
dr 21 

1^2 ll + 

d\ 

dr 22 

\5r 2 2\ 




dO 

1 






dr 23 

y<- 

''23 


dtp 


»3 3 

8 <f> 


r l3 . 

dr i3 


i - 

-'‘23 

: 

= r 

~rW 

d<P 


f 

- 22 

dip 


r 2 l 

dr 2 1 



-'23 

' dr 22 

= T 

— r 23 


(3.7a) 


(3.7b) 


(3.7c) 


(3.7d) 


(3.7e) 

(3.70 


IV. Experimental Results 

In this section, we present experimental results for the IPC and the GIPC 
algorithm tested successfully on real data. The various error plots will also be 
discussed. These plots indicate the relationship between the errors in the 
input set of data (coordinates of the features) and the errors in the output 
data (motion parameters). 

4.1. Using Real Data 


Separate experiments, corresponding to 2-D images of the three positions 
of an octbox in Fig. 3(a), were conducted. The octbox has two parallel 
octagonal faces opposite to each other, and eight rectangular faces. In the first 
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(c) 


Fig. 3. Experiments with real data, (a) Octbox in its initial position; (b) first case of motion 
where Octbox is rotated through 15° around x-axis; and (c) second case of motion where Octbox 
is rotated through 105° around x-axis. 


experiment (Fig. 3(b)), the octbox was rotated around the x-axis by - 15°. In 
the second experiment (Fig. 3(c)), the octbox was rotated around the x-axis 
by 105°. In these experiments, rotation aroudn the x-axis means that the 
angle of rotation, by definition, is roll , or equivalently, the direction cosines of 
the arbitrary axis, around which the octbox rotates, are given by Vj = 1.0; 
v 2 = 0.0; v 3 = 0.0. The translation along the three axes is 1 unit each. The 
data for these three cases of octbox rotations are shown in Figs. 4{a) and 4(b), 
respectively, where coordinates of the vertices of the octbox before and after 
the motion are given. 


4.1.1 . Data Acquisition and Digitization 

A solid octbox with side dimensions 1.5" x 4.5" was constructed and 
placed on a mount capable of motion in all three axes. The mount was flat 
black to minimize reflection from it A video camera was placed at the same 
height above the floor as the octbox and focused on it. The illumination was 
provided by a television floodlight also at the same height as the octbox. The 
octbos was videotaped in 30-s intervals with the illumination source moved 
from being placed along the axis of the camera to a 45° angle from the camera 
and finally to a right angle from the camera. These scenes were recorded on a 
VHS videotape recorder on standard tape at the fastest tape speed allowed. 
The tape was taken to the Image Processing/Computer Vision Laboratory at 
Rice University for digitizing and processing. The tape was digitized using a 
Chorus Data Systems PC-EYE digitizer. Mounted in an IBM-XT, this 
digitizer is capable of capturing a 640 (H) x 400 (V)pixel image with 6 bits 
per pixel quantization. The images were then transferred to the main image 
processing computer system for enhancement and analysis. 
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(b) 



Fig. 4. Real data for motion of Octbox. (a) Data for the first case of motion, and (b) data for 
the second ease of motion. 


4.1.2. Wireframe Extraction 

Our procedure for wireframe extraction consisted of processing the digi- 
tized picture for noise removal and edge detection. In both the raw pictures, 
noise was removed by Wiener filtering. The transfer function is of the form 

™ ^ /**(«, v) 

^“i (4,) 

In our experiment, we assumed that the degradation process //(u, p) was a 
low-pass filter and the noise-to-signal ratio v)/Sf/u y p) was equal to 
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2a 2 . The procedure for removing the noise consisted of the following steps: 

Step 1. Separate the processed picture into tSvo regions according to its 
intensity level. 

Step 2. Compute the variance of each region. 

Step 3. Process each region using Eq. (4.1) with these variances. 

For edge detection, the Sobel edge detector was used. From edges, 
information about the coordinates of the corners was obtained by comparing 
the magnitudes of the gradients with a preset threshold value. 

The GIPC algorithm has been applied to the first and second experiments, 
and the results are shown in Figs. 5(a) and 5(b). In both the experiments, the 
camera was rotated through 10° around its x-axis. With the same set of data 


Estimated Translational Vector ( up to a scale factor) is: 
T - [ 3.863706. 3.863707, 3.863717 ] T 


Two possible solutions of Rotation Matrix are: 

.333334 0.816496 0.471405] 1.000000 0.000000 -0.000000 

1.772304 -0.050319 0.633257 and R' = -0.000000 0.996 1 95 -0.087 156 

1.540773 0.575154 -0.61381QJ 0.000000 0.087156 0.996195 


The directional cosines of the axis and the angle of rotation about the 
axis (corresponding to R and R 1 ) are respectively: 
vj = 0.576984; v 2 = 0.688845; v 3 - 0.438843; 8 = 182.886128 
vf = 1.000000; v 2 ’ « 0.000006; v 3 ' * 0.000001; 0* * 354.999983 
Conclusion : Choose R' and its associated parameters as the final solution. 


Estimated Translational Vector ( up to a scale factor) is: 
T * [ 1.035277, 1.035273, 1.035283 ] T 


Two possible solutions of Rotation Matrix are: 

1.000000 -0.000002 0.000003] [-0.333332 0.471405 -0.816496 

R= 0.000003 -0.087156 -0.996195 and R'= 0.772304 0.633257 0.050321 

0.000002 0.996195 -0.087I56J 0.540773 -0.613810 -0.575153 


The directional cosines of the axis and the angle of rotation about the 
axis (corresponding to R and R*) are respectively: 
v t = 1.000000; v 2 *0.000000; v 3 * 0.000002; 6 * 95.000000 
v i * * 0.431055; v 2 * * 0.860937; v 3 ' * -0.195299; O' * 230.385837 
Conclusion'. Choose R and its associated parameters as the final solution. 


Fig. 5. Demonstration of GIPC algorithm using real data, (a) For the first case of Octbox 
rotation, and (b) for the second case of Octbox rotation. 
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for these experiments for the octbox rotation (Figs. 3(b), 3(c)), the directional 
cosines are found to be same. The angles of rotation are —5° and 95° 
respectively, which means that the angles of rotation of the octbox are added 
to by the amount of rotation by the camera, an9 that indeed should be the 
case. These two experiments show the success of GIPC algorithm with real 


V. Concluding Remarks 

In this paper, a generalized expression for motion-analysis equation was 
derived, and the other three cases of motion-analysis were found to be special 
cases of this case. In addition, the expression for error bounds were dervied 
for various stages of the algorithm. 
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